{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Chapter 20 Non-linear Function Optimization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "from numpy import array, dot, finfo, log, mean, pi, sqrt, squeeze, var, zeros\n",
    "from numpy.linalg import lstsq\n",
    "from numpy.random import randn\n",
    "from scipy.optimize import fmin, fmin_bfgs, fmin_slsqp\n",
    "\n",
    "\n",
    "# Negative Log-likellihood\n",
    "def normal_loglik(mu, x):\n",
    "    n = x.size\n",
    "    epe = dot(x - mu, x - mu)\n",
    "    sigma2 = epe / n\n",
    "    return 0.5 * (n * log(2 * pi) + n * log(sigma2) + epe / sigma2)\n",
    "\n",
    "\n",
    "# Negative Log-likellihood 2\n",
    "def normal_loglik2(p, x):\n",
    "    mu = p[0]\n",
    "    sigma2 = p[1]\n",
    "    if sigma2 < 0:\n",
    "        sigma2 = finfo(float).tiny\n",
    "    epe = dot(x - mu, x - mu)\n",
    "    n = x.size\n",
    "    return 0.5 * (n * log(2 * pi) + n * log(sigma2) + epe / sigma2)\n",
    "\n",
    "\n",
    "# Negative Log-likellihood 3\n",
    "def normal_loglik3(p, x):\n",
    "    mu = p[0]\n",
    "    sigma2 = p[1] ** 2\n",
    "    epe = dot(x - mu, x - mu)\n",
    "    n = x.size\n",
    "    return 0.5 * (n * log(2 * pi) + n * log(sigma2) + epe / sigma2)\n",
    "\n",
    "\n",
    "# Negative Log-likellihood 4\n",
    "def ols_loglik(b, y, X):\n",
    "    bv = b.view()\n",
    "    K = X.shape[1]\n",
    "    bv.shape = K, 1\n",
    "    e = y - dot(X, bv)\n",
    "    n = e.size\n",
    "    epe = dot(e.T, e)\n",
    "    sigma2 = epe / n\n",
    "    return squeeze(0.5 * (n * log(2 * pi) + n * log(sigma2) + epe / sigma2))\n",
    "\n",
    "\n",
    "# Constraint function\n",
    "def constraint(p, x):\n",
    "    return array([1, p[1] - 0.0001])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 1    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = randn(1000) + 2\n",
    "mu0 = mean(x) + 1\n",
    "args = (x,)\n",
    "mu = fmin_bfgs(normal_loglik, mu0, args=args)\n",
    "print(\"mu:\")\n",
    "print(mu)\n",
    "print(\"mean(x):\")\n",
    "print(mean(x))\n",
    "mu = fmin(normal_loglik, mu0, args=args)\n",
    "print(\"mu:\")\n",
    "print(mu)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 2\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p0 = zeros(2)\n",
    "p0[0] = mean(x) + 1\n",
    "p0[1] = var(x) + 2\n",
    "print(\"normal_loglik2(p0,x):\")\n",
    "print(normal_loglik2(p0, x))\n",
    "p = fmin_slsqp(normal_loglik2, p0, args=args, bounds=[(-100, 100), (0.000001, 100)])\n",
    "print(\"array([mean(x),var(x)]):\")\n",
    "print(array([mean(x), var(x)]))\n",
    "print(\"p:\")\n",
    "print(p)\n",
    "p = fmin_slsqp(normal_loglik2, p0, args=args, f_ieqcons=constraint)\n",
    "print(\"p:\")\n",
    "print(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## Exercise 3\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "from scipy.optimize import minimize\n",
    "\n",
    "args = (x,)\n",
    "p = minimize(normal_loglik2, p0, args=args, bounds=[(-100, 100), (0.000001, 100)])\n",
    "print(\"array([mean(x),var(x)]):\")\n",
    "print(array([mean(x), var(x)]))\n",
    "print(\"p:\")\n",
    "print(p)\n",
    "minimize_cons = {\"type\": \"ineq\", \"fun\": lambda p: p[1] - 0.0001}\n",
    "p = minimize(normal_loglik2, p0, args=args, constraints=minimize_cons)\n",
    "print(\"p:\")\n",
    "print(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 4\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p0[1] = sqrt(p0[1])\n",
    "p = fmin_bfgs(normal_loglik3, p0, args=args)\n",
    "p[1] = p[1] ** 2\n",
    "print(\"p:\")\n",
    "print(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 5\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b = array([[1.0], [2]])\n",
    "X = randn(1000, 2)\n",
    "y = dot(X, b) + randn(1000, 1)\n",
    "b0 = lstsq(X, y, rcond=None)\n",
    "b0 = b0[0]\n",
    "print(\"b0:\")\n",
    "print(b0)\n",
    "args = (y, X)\n",
    "print(\"ols_loglik(b0.ravel(),y,X):\")\n",
    "print(ols_loglik(b0.ravel(), y, X))\n",
    "b = fmin_slsqp(ols_loglik, b0.ravel(), args=args)\n",
    "print(\"b:\")\n",
    "print(b)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
